This file analyzes the changes in S/H cell ratio (bleaching) in Pocillopora sp: * In 29 colonies sampled in 2014 (Before ENSO), 2015 (during 1st peak of ENSO) and 2016 (during 2nd peak of ENSO) * These samples were collected at Uva Island, Gulf of Chiriquí, Panama

# Libraries    
    library(tidyverse) # install.packages('tidyverse')
    # library(devtools)  # install.packages("devtools")
    # devtools::install_github("jrcunning/steponeR")
    library(steponeR)
    library(scales)
    library(gridExtra)
    library(lme4)
    library(lmerTest)
    library(multcomp)
    library(emmeans)
   
# Default ggplot settings
    ggthe_bw<-theme(plot.background=element_blank(),
          panel.grid.major.y = element_blank(),
          panel.grid.major.x = element_blank(),
          panel.grid.minor.x = element_blank(),
          panel.grid.minor.y = element_blank(),
          legend.box.background = element_rect(),
          panel.background =element_rect(fill = NA, color = "black")
          )+
    theme_bw()

1. Symbiont to host (S/H) cell ratio and Durusdinium proportion

Get the raw data for all Pocillopora samples applying R.Cunning steponeR function

  1. Get list of plate files to read data
  2. Calculate the ratios
  3. Export results

Get colony and year information from sample labels

# Parse sample names, times, treatments, colonies, plates
    sample.year <- rbind.fill(lapply(strsplit(as.character(Pdam$Sample.Name), split="_"), 
     function(X) data.frame(t(X))))
        colnames(sample.year) <- c("Spp", "Colony","Year")
        Pdam <- cbind(sample.year, Pdam[,-1])

# Keep unique sample ID to check reruns  
    Pdam$Sample<-paste(Pdam$Year,Pdam$Colony, sep="_")
    Pdam$ID<-paste(Pdam$Sample,Pdam$File.Name, sep=".")

Calculate total S/H ratio and D/C ratio and propD

Most of the colonies are dominated by one symbiont genus

2. qPCR DATA CLEANING and quality control

 ## 1. Check and remove NTC wells
  ntc <- Pdam[which(Pdam$Colony=="NTC"), ]
  Pdam <- droplevels(Pdam[!rownames(Pdam) %in% rownames(ntc), ])
  
  ## 2. Check and remove NoHSratio samples
  NoHSratio <- Pdam[which(Pdam$tot.SH==0), ]
  Pdam <- droplevels(Pdam[!rownames(Pdam) %in% rownames(NoHSratio), ])

  ## 3. Chose bw samples ran more than once
  ReRunA <- Pdam[duplicated(Pdam$Sample),]
  
  n_RunA <- data.frame(table(Pdam$Sample))
  colnames(n_RunA)<-c("Sample","RanA")
  Pdam<-join(Pdam, n_RunA, type = "left")
  
  DuplicatesA <- Pdam[(Pdam$RanA>1),]
  
  # Remove bad replicates
  ToRem1<-read.csv("ToRemove10-28-16.csv")  # 11/1/16
  Pdam<-Pdam[!(Pdam$ID %in% ToRem1$ID),]
    
  n_RunB <- data.frame(table(Pdam$Sample))
  ReRunB <- Pdam[duplicated(Pdam$Sample),]
  colnames(n_RunB)<-c("Sample","RanB")
  Pdam<-join(Pdam, n_RunB, type = "left")
  
  # List of dupplicated samples, should have 0 rows now
  DuplicatesB <- Pdam[(Pdam$RanB>1),]
  
  # 4. Check and remove samples with high SD or late amplification
  StDe1.5 <- Pdam[which(Pdam$Pdam.CT.sd>1.5), ]
  Pdam<-Pdam[!(Pdam$ID %in% StDe1.5$ID),]
# End of qPCR data cleaning

3. Field/sample information

Get collection information and select only samples from Uva Island that were sampled in 2014, 2015 and 2016

# 1. Labeling and Factors

  ## Import treatment info from Field file 
  Info<-read.csv("SampleInfo.csv", header=TRUE)

  ## Merge SH ratios and sample info
      Pdam<-join(Pdam, Info, by = "Sample", 
                type = "left", match = "all")
      summary(Pdam)
##      Colony    Year          C.SH               D.SH         
##  256-1  :  3   14:142   Min.   :0.000000   Min.   :0.000000  
##  257-1  :  3   15: 86   1st Qu.:0.000000   1st Qu.:0.000000  
##  259-1  :  3   16: 56   Median :0.006687   Median :0.003965  
##  260-1  :  3            Mean   :0.070606   Mean   :0.034897  
##  262-1  :  3            3rd Qu.:0.036245   3rd Qu.:0.035000  
##  264-1  :  3            Max.   :2.046738   Max.   :1.905453  
##  (Other):266                                                 
##     Sample               ID                tot.SH        
##  Length:284         Length:284         Min.   :0.001125  
##  Class :character   Class :character   1st Qu.:0.018868  
##  Mode  :character   Mode  :character   Median :0.035563  
##                                        Mean   :0.105503  
##                                        3rd Qu.:0.078992  
##                                        Max.   :2.046738  
##                                                          
##    logTot.SH          logC.SH           logD.SH           D.Prp       
##  Min.   :-2.9487   Min.   :-5.5123   Min.   :-5.726   Min.   :0.0000  
##  1st Qu.:-1.7243   1st Qu.:-2.4703   1st Qu.:-1.862   1st Qu.:0.0000  
##  Median :-1.4490   Median :-1.7240   Median :-1.544   Median :0.2061  
##  Mean   :-1.3883   Mean   :-1.9408   Mean   :-1.871   Mean   :0.4674  
##  3rd Qu.:-1.1024   3rd Qu.:-1.2207   3rd Qu.:-1.250   3rd Qu.:1.0000  
##  Max.   : 0.3111   Max.   : 0.3111   Max.   : 0.280   Max.   :1.0000  
##                    NA's   :77        NA's   :114                      
##       RanA            RanB            Region                Location  
##  Min.   :1.000   Min.   :1   Gal_Darwin  : 58   Canal de afuera : 15  
##  1st Qu.:1.000   1st Qu.:1   Gal_Floreana: 12   Contadora Island: 36  
##  Median :1.000   Median :1   Pan_Chiriqui:158   Darwin          : 58  
##  Mean   :1.092   Mean   :1   Pan_Panama  : 56   Floreana        : 12  
##  3rd Qu.:1.000   3rd Qu.:1                      Saboga Island   : 20  
##  Max.   :3.000   Max.   :1                      Secas Island    : 16  
##                                                 Uva Island      :127  
##                          Reef         Depht          Spp     
##  Coral Collection Point    : 40   Min.   : 3.00   P.dam:164  
##  By the Pier               : 34   1st Qu.:10.00   P.eff: 15  
##  Darwin Reef, next to BEAMS: 33   Median :20.00   P.ele: 61  
##  Shallow Millepora         : 24   Mean   :23.24   P.eyd: 14  
##  M.intrincata Shallow      : 17   3rd Qu.:36.00   P.mea: 12  
##  Ridley Rock Shallow       : 16   Max.   :68.00   P.ver: 18  
##  (Other)                   :120                              
##      Year.2         Number        Tag       Collection.day
##  Min.   :2014   Min.   :1    no tag :  9           : 70   
##  1st Qu.:2014   1st Qu.:1    270-1  :  4   08/19/14: 21   
##  Median :2014   Median :1    336    :  4   08/18/15: 17   
##  Mean   :2015   Mean   :1    339-1  :  4   08/27/15: 17   
##  3rd Qu.:2015   3rd Qu.:1    340-1  :  4   04/18/16: 16   
##  Max.   :2016   Max.   :1    345    :  4   04/16/16: 14   
##                 NA's   :78   (Other):255   (Other) :129   
##                                  Notes        Depth.2   
##                                     :274   mid    : 86  
##  1st Rock at 60ft, not the main rock:  3   Shallow:198  
##  No tag                             :  1                
##  No Tag                             :  3                
##  Not 259                            :  1                
##  P.cap                              :  1                
##  Sampled 16 and 17                  :  1
# 2. Uva Data
      # Write how many years a colony was sampled
      Years <- data.frame(table(Pdam$Colony))
      colnames(Years)<-c("Colony","Times")
      Pdam<-join(Pdam, Years, type = "left")
  # Select data only from Uva Island    
      data.Uva<-Pdam[which(Pdam$Location=="Uva Island"), ]
  # Select colonies that were sampled three times (2014, 2015 and 2016)
      data.Uva<-data.Uva[which(data.Uva$Times>2),]
  # Remove "Ross" samples (colony in deep location (~50feet), not Uva reef) 
  # and 259 as well as 264 that were mislabeled in at least one of the sampling points   
      data.Uva<-filter(data.Uva, Colony!="Ross")
      data.Uva<-filter(data.Uva, Colony!="259-1")
      data.Uva<-filter(data.Uva, Colony!="264-1")

      data.Uva<-data.Uva[,c("Colony","Year.2","C.SH","D.SH","tot.SH",
                            "logTot.SH","D.Prp", "logC.SH", "logD.SH","Year","Sample")]
      data.Uva$Year<-factor(data.Uva$Year, levels=c("14", "15", "16"))
      data.Uva$Colony<-factor(data.Uva$Colony)
    summary(data.Uva)
##      Colony       Year.2          C.SH               D.SH        
##  256-1  : 3   Min.   :2014   Min.   :0.000000   Min.   :0.00000  
##  257-1  : 3   1st Qu.:2014   1st Qu.:0.000000   1st Qu.:0.00000  
##  260-1  : 3   Median :2015   Median :0.003005   Median :0.01442  
##  262-1  : 3   Mean   :2015   Mean   :0.019174   Mean   :0.02071  
##  278-1  : 3   3rd Qu.:2016   3rd Qu.:0.011841   3rd Qu.:0.03001  
##  284-1  : 3   Max.   :2016   Max.   :0.690021   Max.   :0.18875  
##  (Other):69                                                      
##      tot.SH          logTot.SH           D.Prp           logC.SH       
##  Min.   :0.00146   Min.   :-2.8356   Min.   :0.0000   Min.   :-5.0627  
##  1st Qu.:0.01294   1st Qu.:-1.8888   1st Qu.:0.0000   1st Qu.:-2.5001  
##  Median :0.02266   Median :-1.6448   Median :0.8620   Median :-2.0971  
##  Mean   :0.03988   Mean   :-1.6539   Mean   :0.5411   Mean   :-2.1808  
##  3rd Qu.:0.03579   3rd Qu.:-1.4463   3rd Qu.:1.0000   3rd Qu.:-1.7015  
##  Max.   :0.69002   Max.   :-0.1611   Max.   :1.0000   Max.   :-0.1611  
##                                                       NA's   :32       
##     logD.SH        Year       Sample         
##  Min.   :-5.7262   14:29   Length:87         
##  1st Qu.:-1.8381   15:29   Class :character  
##  Median :-1.6245   16:29   Mode  :character  
##  Mean   :-2.0059                             
##  3rd Qu.:-1.4640                             
##  Max.   :-0.7241                             
##  NA's   :29

4. Community classification by genera presence/dominance

Symbiodiniaceae genera detected in the samples (Only-C, Only-D, Mixed communities)

# Add a column for the genera detected using qPCR
  data.Uva$Community<-NULL
  data.Uva$Community[which(data.Uva$D.Prp>=0)] <- "Mixed"
  data.Uva$Community[which(data.Uva$D.Prp==0)] <- "Cladocopium-only"
  data.Uva$Community[which(data.Uva$D.Prp==1)] <- "Durusdinium-only"

  data.Uva$Community<-factor(as.character(data.Uva$Community), 
                         levels=c("Cladocopium-only", "Mixed", "Durusdinium-only"))
  
# Summary genera detected
  Summary_Community<-data.Uva %>%
        group_by(Year.2, Community) %>%
        dplyr::summarise(total.count=n()) 
  Summary_Community

Dominant Symbiodiniaceae genus in each sample (C-dominated, D-dominated)

In general, Pocillipora colonies are dominated by one Symbiodiniaceae type that accounts for more than 90 of the symbiont community. However, sometimes intermediate proportions are common when the symbiont community is transitioning from being dominated by one algal type to other (e.g. during heat stress).

For this data set, 2 of the colonies sampled in 2014 had relativelly even proportions of Cladocopium_b and Durusdinim glynnii (Colony 267_2014, D proportion = 0.47; Colony 269_2014, D proportion = 0.50). This made difficult their classification into one domminant type for statistical testing and further plotting. Since these 2 colonies were later dominated by Durusdinim glynnii in 2015, we decided to include them in the D-dominated classification, eventhoug in 2014 there were not such a dominance (See figure 2).

 # Add a column for the dominant genus
  data.Uva$DominantCommunity<-NULL
  #data.Uva$Community[which(data.Uva$D.Prp>0.1 & data.Uva$D.Prp<0.50)] <- "CD"
  #data.Uva$Community[which(data.Uva$D.Prp>0.5 & data.Uva$D.Prp<0.9)] <- "DC"
  data.Uva$DominantCommunity[which(data.Uva$D.Prp>=0.45)] <- "Durusdinium-dominated"  
  data.Uva$DominantCommunity[which(data.Uva$D.Prp<0.45)] <- "Cladocopium-dominated"
  
  data.Uva$DominantCommunity<-factor(as.character(data.Uva$DominantCommunity),
                              levels=c("Cladocopium-dominated","Durusdinium-dominated"))
# Summary dominant genera
  Summary_Community2<-data.Uva %>%
        group_by(Year.2, DominantCommunity) %>%
        dplyr::summarise(total.count=n()) 
  Summary_Community2

Changes in the number of colonies dominated by each genus

To test the statistical significance of changes in the number of colonies dominated by Durusdinium versus Cladocopium we considered inclussion and exclussion of these 2 initially even colonies, and to place them in the D-dominated or C-dominated categories in 2014.

  • Fisher test including initially (Aug 2014) even colonies. Both colonies placed in the D-dominated category.
Input =("
Year    C     D
2014    15   14
2016    8    21
")

Matriz <- as.matrix(read.table(textConnection(Input),
                               header=TRUE, 
                               row.names=1))
fisher.test(Matriz,
            alternative="two.sided")
## 
##  Fisher's Exact Test for Count Data
## 
## data:  Matriz
## p-value = 0.1064
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.8336155 9.7658946
## sample estimates:
## odds ratio 
##   2.760908
  • Fisher test including initially (Aug 2014) even colonies. One colony was placed in the C-dominated category (D proportion=0.47), while the second one was placed in the D-dominated category (D proportion = 0.50)
Input =("
Year    C     D
2014    16   13
2016    8    21
")

Matriz <- as.matrix(read.table(textConnection(Input),
                               header=TRUE, 
                               row.names=1))
fisher.test(Matriz,
            alternative="two.sided")
## 
##  Fisher's Exact Test for Count Data
## 
## data:  Matriz
## p-value = 0.06099
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   0.9556883 11.2396870
## sample estimates:
## odds ratio 
##   3.162767
  • Fisher test removing initially (Aug 2014) even colonies. N=27 colonies instead 29
Input =("
Year    C     D
2014    15   12
2016    8    19
")

Matriz <- as.matrix(read.table(textConnection(Input),
                               header=TRUE, 
                               row.names=1))
fisher.test(Matriz,
            alternative="two.sided")
## 
##  Fisher's Exact Test for Count Data
## 
## data:  Matriz
## p-value = 0.09778
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   0.848568 10.659493
## sample estimates:
## odds ratio 
##   2.906792
  • Fisher test including initially (Aug 2014) even colonies. Both colonies placed in the C-dominated category (D proportion=0.47, colony 267; and 0.50, colony 269)
Input =("
Year    C     D
2014    17   12
2016    8    21
")

Matriz <- as.matrix(read.table(textConnection(Input),
                               header=TRUE, 
                               row.names=1))
fisher.test(Matriz,
            alternative="two.sided")
## 
##  Fisher's Exact Test for Count Data
## 
## data:  Matriz
## p-value = 0.03295
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   1.094473 13.006776
## sample estimates:
## odds ratio 
##   3.629777

5a. qPCR based analysis (Symbiodiniaceae genera hosted)

Set initial (2014) and final (2016) communities for each sample

  # Subset Aug 2014,  Aug 2015 and April 2016 samples
  D2014.data<-subset(data.Uva, Year.2==2014)
  D2014.data<-D2014.data[,c("Colony","D.Prp","Community","DominantCommunity")]
  colnames(D2014.data)<-(c("Colony","D14.D.Prp","D14Community", "D14Dominant"))
  data.Uva<-left_join(data.Uva, D2014.data, by=c("Colony")) # Add initial community to all samples
   
  # D2015.data<-subset(data.Uva, Year.2==2015)
  # D2015.data<-D2015.data[, c("Colony","D.Prp","Community")]
  # colnames(D2015.data)<-(c("Colony","D15.D.Prp","D15Community"))
  # # data.Uva<-left_join(data.Uva, D2015.data, by=c("Colony"))

  D2016.data<-subset(data.Uva, Year.2==2016) # D1a under control temperature
  D2016.data<-D2016.data[, c("Colony","D.Prp","Community","DominantCommunity")]
  colnames(D2016.data)<-(c("Colony","D16.D.Prp","D16Community", "D16Dominant"))
  data.Uva<-left_join(data.Uva, D2016.data, by=c("Colony")) # Add final community to all samples

6a. Changes in colony proportions based on detected genera (qPCR)

Histo_Dprop <- ggplot(data.Uva, aes (D.Prp , fill=(D14Community))) + ggthe_bw +
  geom_histogram(binwidth = 0.05, boundary=0) +
  facet_grid(~Year.2, labeller = labeller(Year.2=c(
                `2014` = "August 2014", `2015` = "August 2015", `2016` = "April 2016"))) +
  scale_fill_manual(values=c("blue3","purple", "red3"),
                    guide = guide_legend(title.position = "top", nrow = 1),
                    name = "Initial symbiont community (August 2014)", 
                    labels=c("Cladocopium-only", "Mixed", "Durusdinium-only")) +
    scale_x_continuous(name= expression(
    italic(Durusdinium)~proportion~(D/H)/(S/H)),
                     breaks=(seq(0,1,by=0.1))) +
    scale_y_continuous((name= "Number of colonies"),
                     breaks=seq(0,15,by=5)) +
    theme(legend.position="bottom",
        legend.title.align=0.5,
        strip.background =element_rect(fill=NA))
Histo_Dprop +coord_flip()

Traject_Proportion_Facet <- ggplot (data.Uva) +
   ggthe_bw + 
    geom_line(aes (x=Year.2, D.Prp, colour=factor(Colony)))+ 
    facet_wrap(~D14Community, labeller = labeller(
    D14Community=c(`Cladocopium-only` = "Cladocopium-only (2014)", 
                   `Mixed` = "Mixed (2014)", 
                   `Durusdinium-only` = "Durusdinium-only (2014)"))) +
    geom_jitter(aes(Year.2, D.Prp, colour=factor(Colony)), alpha=0.3, 
              height = 0.05, width = 0.15) +
    #labs(x = "Año del muestreo", y = "Durusdinium Proportion") +
    theme(legend.position="none")
Traject_Proportion_Facet

Figure 2: qPCR

Histo2014_Dprop <- ggplot(data.Uva[which(data.Uva$Year.2==2014), ], aes (D.Prp , fill=(D14Community))) +
      ggthe_bw + coord_flip()+ ggtitle("a.")+
      geom_histogram(binwidth = 0.05, boundary=0) +
      scale_fill_manual(values=c("blue3","purple", "red3"),
                        guide = guide_legend(title.position = "top", nrow = 1),
                        name = "Initial symbiont community (August 2014)", 
                        labels=c("Cladocopium-only", "Mixed", "Durusdinium-only")) +
      scale_x_continuous(name= expression(
                        italic(Durusdinium)~proportion~((D/H)/(S/H))),
                        breaks=(seq(0,1,by=0.1)),
                        expand=c(0.01,0.01)) +
      scale_y_continuous((name= "Number of colonies \n (Aug 2014)"),
                        limits = c(0,17),
                        breaks=seq(0,15,by=5)) +
      theme(legend.position="bottom",
            legend.title.align=0.5,
            strip.background =element_rect(fill=NA))
#Histo2014_Dprop 

Histo2016_Dprop <- ggplot(data.Uva[which(data.Uva$Year.2==2016), ], aes (D.Prp , fill=(D14Community))) +
      ggthe_bw + coord_flip()+ ggtitle("c.")+
      geom_histogram(binwidth = 0.05, boundary=0) +
      scale_fill_manual(values=c("blue3","purple", "red3"),
                        guide = guide_legend(title.position = "top", nrow = 1),
                        name = "Initial symbiont community (August 2014)", 
                        labels=c("Cladocopium-only", "Mixed", "Durusdinium-only")) +
      scale_x_continuous(name="",
                        breaks=(seq(0,1,by=0.1)),
                        expand=c(0.01,0.01)) +    
      scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
                        limits = c(0,17),
                        breaks=seq(0,15,by=5)) +
      theme(legend.position="bottom",
            legend.title.align=0.5,
            strip.background =element_rect(fill=NA))
#Histo2016_Dprop

Traject_Proportion<- ggplot(data.Uva, aes (x=Year.2, y=D.Prp , colour=(D14Community), group=(Colony))) + ggthe_bw + ggtitle("b.") +
  geom_line(linetype=2) +
  geom_jitter(height = 0.01, width = 0.25, alpha=0.4) +
  
  scale_colour_manual(values=c("blue3","purple", "red3"),
                    guide = guide_legend(title.position = "top", nrow = 1),
                    name = "Initial symbiont community (August 2014)", 
                    labels=c("Cladocopium-only", "Mixed", "Durusdinium-only")) +
  scale_x_continuous(name="", 
                         breaks = c(2014, 2015, 2016),
                         labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) +
  scale_y_continuous((name= ""),
                     breaks=seq(0,1,by=0.1),
                     expand=c(0.01,0.01)) +
  
  theme(legend.position="bottom",
        legend.title.align=0.5,
        strip.background =element_rect(fill=NA))
#Traject_Proportion

Figure_2<-grid.arrange(Histo2014_Dprop,Traject_Proportion, Histo2016_Dprop,
nrow=1, widths=c(3/9, 3/9, 3/9))

# ggsave(file="Outputs/Figure_2qPCR.svg", plot=Figure_2, width=7, height=4)

7a. Changes in total Symbiodiniaceae abundance (bleaching): qPCR

SH_2014Co<-lmerTest::lmer(logTot.SH ~ Year * D14Community + (1|Colony), data=data.Uva)
    summary(SH_2014Co)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * D14Community + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 83.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6538 -0.5202 -0.1347  0.3700  2.9430 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.03036  0.1742  
##  Residual             0.10764  0.3281  
## Number of obs: 87, groups:  Colony, 29
## 
## Fixed effects:
##                                     Estimate Std. Error      df t value
## (Intercept)                          -1.6402     0.1238 71.1153 -13.246
## Year15                               -0.4132     0.1547 52.0000  -2.672
## Year16                               -0.5536     0.1547 52.0000  -3.579
## D14CommunityMixed                    -0.0137     0.1638 71.1153  -0.084
## D14CommunityDurusdinium-only          0.1232     0.1805 71.1153   0.682
## Year15:D14CommunityMixed              0.7731     0.2046 52.0000   3.778
## Year16:D14CommunityMixed              0.6960     0.2046 52.0000   3.402
## Year15:D14CommunityDurusdinium-only   0.4584     0.2255 52.0000   2.033
## Year16:D14CommunityDurusdinium-only   0.3869     0.2255 52.0000   1.716
##                                     Pr(>|t|)    
## (Intercept)                          < 2e-16 ***
## Year15                              0.010053 *  
## Year16                              0.000756 ***
## D14CommunityMixed                   0.933599    
## D14CommunityDurusdinium-only        0.497283    
## Year15:D14CommunityMixed            0.000408 ***
## Year16:D14CommunityMixed            0.001294 ** 
## Year15:D14CommunityDurusdinium-only 0.047152 *  
## Year16:D14CommunityDurusdinium-only 0.092125 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Year15 Year16 D14CmM D14CD- Y15:D14CM Y16:D14CM
## Year15      -0.624                                                
## Year16      -0.624  0.500                                         
## D14CmmntyMx -0.756  0.472  0.472                                  
## D14CmmntyD- -0.686  0.428  0.428  0.519                           
## Yr15:D14CmM  0.472 -0.756 -0.378 -0.624 -0.324                    
## Yr16:D14CmM  0.472 -0.378 -0.756 -0.624 -0.324  0.500             
## Yr15:D14CD-  0.428 -0.686 -0.343 -0.324 -0.624  0.519     0.259   
## Yr16:D14CD-  0.428 -0.343 -0.686 -0.324 -0.624  0.259     0.519   
##             Y15:D14CD
## Year15               
## Year16               
## D14CmmntyMx          
## D14CmmntyD-          
## Yr15:D14CmM          
## Yr16:D14CmM          
## Yr15:D14CD-          
## Yr16:D14CD-  0.500
    #step(SH_2014Co)
    
    SH.emmc<-emmeans(SH_2014Co, ~Year*D14Community)
    SH_groups<-multcomp::cld(SH.emmc)
    SH_groups
    plot(emmeans(SH_2014Co, ~Year|~D14Community), comparisons = TRUE) +
      coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_wrap(~D14Community)

SH_2014Com <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(D14Community))) +  ggthe_bw +
      stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1, 
                       position=position_dodge(width=0.95))+
      stat_summary(fun.data = "mean_cl_boot",geom = "point", 
                       position=position_dodge(width=0.95))+
      stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.95)) +
      
      scale_colour_manual(values=c("blue","purple", "red"),
                            guide = guide_legend(title.position = "top", nrow = 1),
                            name = "Symbionts detectec in 2014",
                            labels=c("Cladocopium", "Mixed", "Durusdinium"))+
      theme(legend.position=c(0.5,0.1),
            strip.background =element_rect(fill="white")) +
      scale_y_continuous(limits = c(-3, -0.8),
                name="\n Total symbiont to host cell ratio (S/H)", 
                labels = (math_format(10^.x)),
                expand = c(0.0, 0)) + 
      scale_x_discrete(name="Year")
SH_2014Com

SH_Community<-lmerTest::lmer(logTot.SH ~ Year * Community + (1|Colony), data=data.Uva)
    summary(SH_Community)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * Community + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 91.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3973 -0.5171 -0.0987  0.3905  4.0587 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.02404  0.1550  
##  Residual             0.12526  0.3539  
## Number of obs: 87, groups:  Colony, 29
## 
## Fixed effects:
##                                  Estimate Std. Error       df t value
## (Intercept)                      -1.60009    0.12804 76.71005 -12.497
## Year15                           -0.20778    0.15540 49.73820  -1.337
## Year16                           -0.71036    0.17979 47.25496  -3.951
## CommunityMixed                   -0.09453    0.16829 77.96039  -0.562
## CommunityDurusdinium-only         0.09882    0.18645 77.17623   0.530
## Year15:CommunityMixed             0.58170    0.24003 54.72488   2.423
## Year16:CommunityMixed             0.81742    0.24650 54.96059   3.316
## Year15:CommunityDurusdinium-only  0.27313    0.22958 47.93524   1.190
## Year16:CommunityDurusdinium-only  0.60594    0.24059 49.84009   2.519
##                                  Pr(>|t|)    
## (Intercept)                       < 2e-16 ***
## Year15                           0.187280    
## Year16                           0.000258 ***
## CommunityMixed                   0.575911    
## CommunityDurusdinium-only        0.597640    
## Year15:CommunityMixed            0.018708 *  
## Year16:CommunityMixed            0.001622 ** 
## Year15:CommunityDurusdinium-only 0.240035    
## Year16:CommunityDurusdinium-only 0.015039 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Year15 Year16 CmmntM CmmnD- Y15:CM Y16:CM Y15:CD
## Year15      -0.733                                                 
## Year16      -0.602  0.497                                          
## CommuntyMxd -0.754  0.576  0.457                                   
## CmmntyDrsd- -0.686  0.504  0.414  0.524                            
## Yr15:CmmntM  0.472 -0.658 -0.322 -0.628 -0.326                     
## Yr16:CmmntM  0.459 -0.365 -0.740 -0.608 -0.318  0.412              
## Yr15:CmmnD-  0.496 -0.677 -0.337 -0.383 -0.717  0.440  0.243       
## Yr16:CmmnD-  0.451 -0.369 -0.747 -0.326 -0.682  0.240  0.543  0.541
    #step(SH_Community)
    
    SH.emmc<-emmeans(SH_Community, ~Year*Community)
    SH_groups<-multcomp::cld(SH.emmc)    
    #write.csv(SH_groups,"SH_groupsB.csv")
    
plot(emmeans(SH_Community, ~Year|~Community), comparisons = TRUE) +
      coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_wrap(~Community)

SH_Com <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(Community))) +  ggthe_bw +
      stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1, 
                       position=position_dodge(width=0.95))+
      stat_summary(fun.data = "mean_cl_boot",geom = "point", 
                       position=position_dodge(width=0.95))+
      stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.95)) +
      
      scale_colour_manual(values=c("blue","purple", "red"),
                            guide = guide_legend(title.position = "top", nrow = 1),
                            name = "Symbionts detectec in 2014",
                            labels=c("Cladocopium", "Mixed", "Durusdinium"))+
      theme(legend.position=c(0.5,0.1),
            strip.background =element_rect(fill="white")) +
      scale_y_continuous(limits = c(-3, -0.8),
                name="\n Total symbiont to host cell ratio (S/H)", 
                labels = (math_format(10^.x)),
                expand = c(0.0, 0)) + 
      scale_x_discrete(name="Year")
SH_Com

anova(SH_2014Co, SH_Community)

Figure 3: qPCR

SH_Communities<-lmerTest::lmer(logTot.SH ~ Year * D14Community * Community + (1|Colony), data=data.Uva)
    summary(SH_Communities)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * D14Community * Community + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 81.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.73407 -0.48463 -0.06579  0.35339  2.91935 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.02179  0.1476  
##  Residual             0.11026  0.3320  
## Number of obs: 87, groups:  Colony, 29
## 
## Fixed effects:
##                                     Estimate Std. Error      df t value
## (Intercept)                          -1.6402     0.1211 70.1550 -13.541
## Year15                               -0.4132     0.1565 47.2363  -2.640
## Year16                               -0.7093     0.1688 50.6884  -4.202
## D14CommunityMixed                    -0.7146     0.3265 73.0936  -2.189
## D14CommunityDurusdinium-only         -0.6296     0.3928 72.9975  -1.603
## CommunityMixed                        0.7009     0.2845 70.5547   2.464
## CommunityDurusdinium-only             0.7527     0.3508 70.9827   2.146
## Year15:D14CommunityMixed              1.5019     0.3806 64.8657   3.946
## Year16:D14CommunityMixed              0.8259     0.2396 53.9251   3.447
## Year15:D14CommunityDurusdinium-only   1.2674     0.5110 61.4046   2.480
## Year16:D14CommunityDurusdinium-only   0.5426     0.2368 48.9752   2.292
## Year15:CommunityMixed                -0.7379     0.3653 70.6556  -2.020
## Year15:CommunityDurusdinium-only     -0.8090     0.4572 64.7392  -1.770
##                                     Pr(>|t|)    
## (Intercept)                          < 2e-16 ***
## Year15                              0.011206 *  
## Year16                              0.000107 ***
## D14CommunityMixed                   0.031805 *  
## D14CommunityDurusdinium-only        0.113266    
## CommunityMixed                      0.016178 *  
## CommunityDurusdinium-only           0.035332 *  
## Year15:D14CommunityMixed            0.000198 ***
## Year16:D14CommunityMixed            0.001108 ** 
## Year15:D14CommunityDurusdinium-only 0.015876 *  
## Year16:D14CommunityDurusdinium-only 0.026258 *  
## Year15:CommunityMixed               0.047162 *  
## Year15:CommunityDurusdinium-only    0.081500 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 14 columns / coefficients
    #step(SH_Communities)
    
    SH.emmc<-emmeans(SH_Communities, ~Year*D14Community*Community)
    SH_groups<-cld(SH.emmc)
    SH_groups<-as.data.frame(SH_groups[complete.cases(SH_groups),])
    SH_groups<-SH_groups[order(SH_groups$Year,SH_groups$D14Community, SH_groups$Community),]
    SH_groups
    #write.csv(SH_groups,"SH_Communities_multicomp.csv")

SH_Comties <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(Community))) + ggthe_bw+
    
      facet_grid(~D14Community, labeller = labeller(D14Community=c(
                      `Cladocopium-only` = "Cladocopium-only (2014)",
                      `Mixed` = "Mixed (2014)", 
                      `Durusdinium-only` = "Durusdinium-only (2014)"))) +
        
      stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
      stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
      stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
      
      scale_colour_manual(values=c("blue3","purple", "red3"),
                        guide = guide_legend(title.position = "top", nrow = 3),
                        name = "Detected symbionts",
                        labels=c("Cladocopium", "Mixed", "Durusdinium"))+
      theme(legend.position=c(0.82,0.2),
            legend.box.background = element_rect(),
            strip.background =element_rect(fill="white"))+
   
      scale_y_continuous(limits = c(-3, -0.8),
            name="\n Total symbiont to host cell ratio (S/H)", 
            labels = (math_format(10^.x)),
            expand = c(0.0, 0)) + 
      scale_x_continuous(name="", 
                         breaks = c(14, 15, 16),
                         labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) 
SH_Comties

#ggsave(file="Figure_3.svg", plot=SH_Comties, width=6, height=4)

8a. Changes in Cladocopium abundance

CH_Com <- ggplot(data.Uva, aes (Year.2, logC.SH, colour=factor(Community))) + ggthe_bw+ ggtitle("a.")+
    
      stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
      stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
      stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
      
      scale_colour_manual(values=c("blue3","purple", "red3"),
                        guide = guide_legend(title.position = "top", nrow = 3),
                        name = "Detected symbionts",
                        labels=c("Cladocopium", "Mixed", "Durusdinium"))+
      theme(legend.position=c(0.82,0.8),
            legend.box.background = element_rect(),
            strip.background =element_rect(fill="white"))+
   
      scale_y_continuous(limits = c(-3, -0.7),
            name="\n Cladocopium to host cell ratio (C/H)", 
            labels = (math_format(10^.x)),
            expand = c(0.0, 0)) + 
      scale_x_continuous(name="", 
                         breaks = c(2014, 2015, 2016),
                         labels = c(" ", " ", " ")) 
#CH_Com

CH_Com2<-CH_Com + facet_grid(~D14Dominant, labeller = labeller(D14Dominant=c(
                   `Cladocopium-dominated` = "Cladocopium- \n dominated (2014)",
                   `Durusdinium-dominated` = "Durusdinium- \n dominated (2014)")))


CH<-lmer(logC.SH ~ Year* Community* D14Dominant + (1|Colony), data=data.Uva)
    summary(CH)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logC.SH ~ Year * Community * D14Dominant + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 113.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.70311 -0.36865 -0.09556  0.50582  1.89250 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.1182   0.3438  
##  Residual             0.4032   0.6350  
## Number of obs: 55, groups:  Colony, 21
## 
## Fixed effects:
##                                         Estimate Std. Error       df
## (Intercept)                             -1.54068    0.23810 44.78323
## Year15                                  -0.25625    0.28031 26.58185
## Year16                                  -0.73614    0.32339 24.71137
## CommunityMixed                          -0.32887    0.37022 44.99572
## D14DominantDurusdinium-dominated        -1.09559    0.41353 43.96389
## Year15:CommunityMixed                    0.59388    0.61658 32.70059
## Year16:CommunityMixed                   -0.03386    0.50241 28.81910
## Year15:D14DominantDurusdinium-dominated -0.54917    0.68606 30.91184
## Year16:D14DominantDurusdinium-dominated  0.88798    0.66184 31.43058
##                                         t value Pr(>|t|)    
## (Intercept)                              -6.471 6.36e-08 ***
## Year15                                   -0.914   0.3688    
## Year16                                   -2.276   0.0318 *  
## CommunityMixed                           -0.888   0.3791    
## D14DominantDurusdinium-dominated         -2.649   0.0112 *  
## Year15:CommunityMixed                     0.963   0.3425    
## Year16:CommunityMixed                    -0.067   0.9467    
## Year15:D14DominantDurusdinium-dominated  -0.800   0.4296    
## Year16:D14DominantDurusdinium-dominated   1.342   0.1893    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Year15 Year16 CmmntM D14DD- Y15:CM Y16:CM Y15:D1
## Year15      -0.715                                                 
## Year16      -0.579  0.492                                          
## CommuntyMxd -0.622  0.502  0.369                                   
## D14DmnntDr- -0.019 -0.038  0.003 -0.537                            
## Yr15:CmmntM  0.318 -0.474 -0.224 -0.512  0.275                     
## Yr16:CmmntM  0.400 -0.331 -0.659 -0.643  0.345  0.357              
## Yr15:D14DD-  0.006  0.017  0.000  0.255 -0.468 -0.705 -0.185       
## Yr16:D14DD- -0.021  0.011  0.011  0.308 -0.509 -0.162 -0.437  0.306
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficients
    #step(CH)

    CH.emmc<-emmeans(CH, ~Year*Community* D14Dominant)
    CH_groups<-cld(CH.emmc)
    CH_groups<-as.data.frame(CH_groups[complete.cases(CH_groups),])
    CH_groups<-CH_groups[order(CH_groups$Year,CH_groups$D14Dominant, CH_groups$Community),]
    CH_groups
    #write.csv(CH_groups,"CH_groups.csv")

9a. Changes in Durusdinium abundance

DH_Com <- ggplot(data.Uva, aes (Year.2, logD.SH, colour=factor(Community))) + 
      ggthe_bw+ ggtitle("b.")+
    
      stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
      stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
      stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
      
      scale_colour_manual(values=c("purple", "red3"),
                        guide = guide_legend(title.position = "top", nrow = 3),
                        name = "Detected symbionts",
                        labels=c("Mixed", "Durusdinium"))+
      theme(legend.position=c(0.82,0.2),
            legend.box.background = element_rect(),
            strip.background =element_rect(fill="white"))+
   
      scale_y_continuous(limits = c(-6, -0.7),
            name="\n Durusdinium to host cell ratio (D/H)", 
            labels = (math_format(10^.x)),
            expand = c(0.0, 0)) + 
      scale_x_continuous(name="", 
                         breaks = c(2014, 2015, 2016),
                         labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016"))
                           
# DH_Com

DH_Com2<-DH_Com + facet_grid(~D14Dominant, labeller = labeller(D14Dominant=c(
                   `Cladocopium-dominated` = "Cladocopium- \n dominated (2014)",
                   `Durusdinium-dominated` = "Durusdinium- \n dominated (2014)")))


DH_LM<-lmer(logD.SH ~ Year* Community* D14Dominant + (1|Colony), data=data.Uva)
    summary(DH_LM)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logD.SH ~ Year * Community * D14Dominant + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 66.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0953 -0.5756  0.0079  0.3844  3.6310 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.0000   0.0000  
##  Residual             0.1687   0.4107  
## Number of obs: 58, groups:  Colony, 22
## 
## Fixed effects:
##                                                            Estimate
## (Intercept)                                                -4.86573
## Year15                                                      0.07049
## Year16                                                      3.21459
## CommunityDurusdinium-only                                   0.27806
## D14DominantDurusdinium-dominated                            3.15614
## Year15:CommunityDurusdinium-only                           -0.35506
## Year16:CommunityDurusdinium-only                           -0.22094
## Year15:D14DominantDurusdinium-dominated                     0.36688
## Year16:D14DominantDurusdinium-dominated                    -3.09697
## CommunityDurusdinium-only:D14DominantDurusdinium-dominated -0.08557
##                                                            Std. Error
## (Intercept)                                                   0.16766
## Year15                                                        0.33532
## Year16                                                        0.23711
## CommunityDurusdinium-only                                     0.50992
## D14DominantDurusdinium-dominated                              0.23711
## Year15:CommunityDurusdinium-only                              0.32897
## Year16:CommunityDurusdinium-only                              0.38416
## Year15:D14DominantDurusdinium-dominated                       0.42745
## Year16:D14DominantDurusdinium-dominated                       0.41068
## CommunityDurusdinium-only:D14DominantDurusdinium-dominated    0.45916
##                                                                  df
## (Intercept)                                                48.00000
## Year15                                                     48.00000
## Year16                                                     48.00000
## CommunityDurusdinium-only                                  48.00000
## D14DominantDurusdinium-dominated                           48.00000
## Year15:CommunityDurusdinium-only                           48.00000
## Year16:CommunityDurusdinium-only                           48.00000
## Year15:D14DominantDurusdinium-dominated                    48.00000
## Year16:D14DominantDurusdinium-dominated                    48.00000
## CommunityDurusdinium-only:D14DominantDurusdinium-dominated 48.00000
##                                                            t value
## (Intercept)                                                -29.021
## Year15                                                       0.210
## Year16                                                      13.557
## CommunityDurusdinium-only                                    0.545
## D14DominantDurusdinium-dominated                            13.311
## Year15:CommunityDurusdinium-only                            -1.079
## Year16:CommunityDurusdinium-only                            -0.575
## Year15:D14DominantDurusdinium-dominated                      0.858
## Year16:D14DominantDurusdinium-dominated                     -7.541
## CommunityDurusdinium-only:D14DominantDurusdinium-dominated  -0.186
##                                                            Pr(>|t|)    
## (Intercept)                                                 < 2e-16 ***
## Year15                                                        0.834    
## Year16                                                      < 2e-16 ***
## CommunityDurusdinium-only                                     0.588    
## D14DominantDurusdinium-dominated                            < 2e-16 ***
## Year15:CommunityDurusdinium-only                              0.286    
## Year16:CommunityDurusdinium-only                              0.568    
## Year15:D14DominantDurusdinium-dominated                       0.395    
## Year16:D14DominantDurusdinium-dominated                    1.09e-09 ***
## CommunityDurusdinium-only:D14DominantDurusdinium-dominated    0.853    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Year15 Year16 CmmnD- D14DD- Y15:CD Y16:CD Y15:D1 Y16:D1
## Year15      -0.500                                                        
## Year16      -0.707  0.354                                                 
## CmmntyDrsd-  0.000  0.000 -0.232                                          
## D14DmnntDr- -0.707  0.354  0.500 -0.232                                   
## Yr15:CmmnD-  0.000  0.000  0.000 -0.293  0.360                            
## Yr16:CmmnD-  0.000  0.000  0.000 -0.753  0.309  0.389                     
## Yr15:D14DD-  0.392 -0.784 -0.277  0.129 -0.555 -0.500 -0.171              
## Yr16:D14DD-  0.408 -0.204 -0.577  0.671 -0.577 -0.208 -0.713  0.320       
## CmD-:D14DD-  0.000  0.000  0.258 -0.900  0.000  0.000  0.558  0.000 -0.596
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
## convergence code: 0
## boundary (singular) fit: see ?isSingular
    #step(DH_LM)

    DH.emmc<-emmeans(DH_LM, ~Year*Community* D14Dominant)
    DH_groups<-cld(DH.emmc)
    DH_groups<-as.data.frame(DH_groups[complete.cases(DH_groups),])
    DH_groups<-DH_groups[order(DH_groups$Year,DH_groups$D14Dominant, DH_groups$Community),]
    DH_groups
    #write.csv(DH_groups,"DH_groups.csv")

Figure 4: qPCR

Figure_4<-grid.arrange(CH_Com2, DH_Com2,nrow=2)

# ggsave(file="Outputs/Figure_4.svg", plot=Figure_4, width=7, height=4)

5b. DGGE based analysis (Symbiodiniaceae ITS2 type hosted)

Community classification by ITS2 dominant profile in 2014

  # Add initial (2014) ITS2 profiles
    ITS2_C_a<-data.frame("Colony"=c("256-1", "257-1", "260-1", "262-1", "330-1", "340-1A"), 
                         "ITS2_14"="C_a")
    # ITS2_C_b<-data.frame("Colony"=c("265-1", "266-1", "267", "268-1", "269-1", "270-1A", "271-1", "285-1", "287-1", "336-1A", "339-1A"), "ITS2_14"="C_b") # 267 and 269 placed under D1
    
    ITS2_C_b<-data.frame("Colony"=c("265-1", "266-1", "268-1",  "270-1A", "271-1", "285-1",
                                    "287-1", "336-1A", "339-1A"), "ITS2_14"="C_b")
    ITS2_C<-rbind(ITS2_C_a, ITS2_C_b)
    
    ITS2_D<-anti_join(D2014.data, ITS2_C, by="Colony")
    ITS2_D$ITS2_14<-"D1"
    ITS2_D<-ITS2_D[, c("Colony", "ITS2_14")]
    
    ITS2<-rbind(ITS2_C, ITS2_D)
    data.Uva<-left_join(data.Uva, ITS2, by=c("Colony")) # Add initial ITS2 profiles to all samples

6b. Changes in colony proportions based on ITS2 profiles (2014) and genus dominance

Histo_Dprop_ITS2 <- ggplot(data.Uva, aes (D.Prp , fill=(ITS2_14)), shape=(D14Community))  + ggthe_bw +
      geom_histogram(binwidth = 0.05, boundary=0) +
      facet_grid(~Year.2, labeller = labeller(Year.2=c(
                    `2014` = "August 2014", `2015` = "August 2015", `2016` = "April 2016"))) +
      scale_fill_manual(values=c("blue3","purple", "red3"),
                        guide = guide_legend(title.position = "top", nrow = 1),
                        name = "Initial ITS2 community (August 2014)", 
                        labels=c("C_a", "C_b", "D1")) +
      scale_x_continuous(name= expression(
                        italic(Durusdinium)~proportion~(D/H)/(S/H)),
                        breaks=(seq(0,1,by=0.1))) +
      scale_y_continuous((name= "Number of colonies"),
                         breaks=seq(0,15,by=5)) +
      theme(legend.position="bottom",
            legend.title.align=0.5,
            strip.background =element_rect(fill=NA))
Histo_Dprop_ITS2 +coord_flip()

Traject_Proportion_FacetITS2 <- ggplot (data.Uva) +
   ggthe_bw + 
    geom_line(aes (x=Year.2, D.Prp, colour=factor(Colony)))+ 
    facet_wrap(~ITS2_14) +
    geom_jitter(aes(Year.2, D.Prp, colour=factor(Colony)), alpha=0.3, 
              height = 0.05, width = 0.15) +
    #labs(x = "Año del muestreo", y = "Durusdinium Proportion") +
    theme(legend.position="none")
Traject_Proportion_FacetITS2

Figure 2: ITS2

Histo2014_DpropITS2 <- ggplot(data.Uva[which(data.Uva$Year.2==2014), ], aes (D.Prp , fill=(ITS2_14))) +
        ggthe_bw + coord_flip()+ ggtitle("a.")+
        geom_histogram(binwidth = 0.05, boundary=0) +
        scale_fill_manual(values=c("blue3","purple", "red3"),
                          guide = guide_legend(title.position = "top", nrow = 1),
                          name = "Initial ITS2 community (August 2014)", 
                          labels=c("C_a", "C_b", "D1")) +
        scale_x_continuous(name= expression(
                          italic(Durusdinium)~proportion~((D/H)/(S/H))),
                          breaks=(seq(0,1,by=0.1)),
                          expand=c(0.01,0.01)) +
        scale_y_continuous((name= "Number of colonies \n (Aug 2014)"),
                          limits = c(0,17),
                          breaks=seq(0,15,by=5)) +
        theme(legend.position="bottom",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA))
#Histo2014_DpropITS2 

Histo2016_DpropITS2 <- ggplot(data.Uva[which(data.Uva$Year.2==2016), ], aes (D.Prp , fill=(ITS2_14))) +
        ggthe_bw + coord_flip()+ ggtitle("c.")+
        geom_histogram(binwidth = 0.05, boundary=0) +
        scale_fill_manual(values=c("blue3","purple", "red3"),
                          guide = guide_legend(title.position = "top", nrow = 1),
                          name = "Initial ITS2 community (August 2014)", 
                          labels=c("C_a", "C_b", "D1")) +
        scale_x_continuous(name="",
                           breaks=(seq(0,1,by=0.1)),
                           expand=c(0.01,0.01)) +    
        scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
                          limits = c(0,17),
                          breaks=seq(0,15,by=5)) +
        theme(legend.position="bottom",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA))
#Histo2016_DpropITS2

Traject_ProportionITS2<- ggplot(data.Uva, aes (x=Year.2, y=D.Prp , colour=(ITS2_14), group=(Colony)))+
        ggthe_bw + ggtitle("b.") +
        geom_line(linetype=2) +
        geom_jitter(height = 0.01, width = 0.25, alpha=0.4) +
        
        scale_colour_manual(values=c("blue3","purple", "red3"),
                          guide = guide_legend(title.position = "top", nrow = 1),
                          name = "Initial ITS2 community (August 2014)", 
                          labels=c("C_a", "C_b", "D1")) +
        scale_x_continuous(name="", 
                           breaks = c(2014, 2015, 2016),
                           labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) +
        scale_y_continuous((name= ""),
                           breaks=seq(0,1,by=0.1),
                           expand=c(0.01,0.01)) +
        
        theme(legend.position="bottom",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA))
#Traject_ProportionITS2


Figure_2_ITS2<-grid.arrange(Histo2014_DpropITS2, Traject_ProportionITS2,Histo2016_DpropITS2,
nrow=1, widths=c(3/9, 3/9, 3/9))

# ggsave(file="Outputs/Figure_ITS2.svg", plot=Figure_2_ITS2, width=7, height=4)

7b. Changes in total Symbiodiniaceae abundance (bleaching): qPCR + Initial ITS2

Only considering initial (2014) community

SH_2014CoI<-lmerTest::lmer(logTot.SH ~ Year * ITS2_14 + (1|Colony), data=data.Uva)
    summary(SH_2014CoI)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * ITS2_14 + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 77.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3417 -0.5054 -0.0500  0.3438  3.4861 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.0190   0.1378  
##  Residual             0.1058   0.3253  
## Number of obs: 87, groups:  Colony, 29
## 
## Fixed effects:
##                   Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)       -1.66014    0.14424 74.54698 -11.510  < 2e-16 ***
## Year15            -0.56741    0.18782 52.00000  -3.021 0.003902 ** 
## Year16            -0.77804    0.18782 52.00000  -4.142 0.000127 ***
## ITS2_14C_b        -0.01746    0.18621 74.54698  -0.094 0.925547    
## ITS2_14D1          0.11106    0.17240 74.54698   0.644 0.521417    
## Year15:ITS2_14C_b  0.79408    0.24248 52.00000   3.275 0.001885 ** 
## Year16:ITS2_14C_b  0.88586    0.24248 52.00000   3.653 0.000602 ***
## Year15:ITS2_14D1   0.73349    0.22449 52.00000   3.267 0.001926 ** 
## Year16:ITS2_14D1   0.71310    0.22449 52.00000   3.177 0.002508 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##              (Intr) Year15 Year16 ITS2_14C ITS2_14D Y15:ITS2_14C
## Year15       -0.651                                             
## Year16       -0.651  0.500                                      
## ITS2_14C_b   -0.775  0.504  0.504                               
## ITS2_14D1    -0.837  0.545  0.545  0.648                        
## Y15:ITS2_14C  0.504 -0.775 -0.387 -0.651   -0.422               
## Y16:ITS2_14C  0.504 -0.387 -0.775 -0.651   -0.422    0.500      
## Y15:ITS2_14D  0.545 -0.837 -0.418 -0.422   -0.651    0.648      
## Y16:ITS2_14D  0.545 -0.418 -0.837 -0.422   -0.651    0.324      
##              Y16:ITS2_14C Y15:ITS2_14D
## Year15                                
## Year16                                
## ITS2_14C_b                            
## ITS2_14D1                             
## Y15:ITS2_14C                          
## Y16:ITS2_14C                          
## Y15:ITS2_14D  0.324                   
## Y16:ITS2_14D  0.648        0.500
    #step(SH_2014CoI)
    
    SH_I.emmc<-emmeans(SH_2014CoI, ~Year*ITS2_14)
    SH_I_groups<-multcomp::cld(SH_I.emmc)
    SH_I_groups
    SH_I_groups<-as.data.frame(SH_I_groups[complete.cases(SH_I_groups),])
    SH_I_groups<-SH_I_groups[order(SH_I_groups$Year,SH_I_groups$ITS2_14),]
    SH_I_groups
    plot(emmeans(SH_2014CoI, ~Year|~ITS2_14), comparisons = TRUE) +
      coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_wrap(~ITS2_14)

SH_2014ComI <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(ITS2_14))) +  ggthe_bw +
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1, 
                   position=position_dodge(width=0.95))+
  stat_summary(fun.data = "mean_cl_boot",geom = "point", 
                   position=position_dodge(width=0.95))+
  stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.95)) +
  
  scale_colour_manual(values=c("blue","purple", "red"),
                        guide = guide_legend(title.position = "top", nrow = 1),
                        name = "ITS2 detected in 2014",
                        labels=c("C_a", "C_b", "D1"))+
  theme(legend.position=c(0.5,0.1),
        strip.background =element_rect(fill="white")) +
  scale_y_continuous(limits = c(-3, -0.8),
            name="\n Total symbiont to host cell ratio (S/H)", 
            labels = (math_format(10^.x)),
            expand = c(0.0, 0)) + 
  scale_x_discrete(name="Year")
SH_2014ComI

Only considering current community ADD CURRENT ITS2!

# SH_Community<-lmerTest::lmer(logTot.SH ~ Year * Community + (1|Colony), data=data.Uva)
#     summary(SH_Community)
#     #step(SH_Community)
#     
#     SH.emmc<-emmeans(SH_Community, ~Year*Community)
#     SH_groups<-multcomp::cld(SH.emmc)    
#     #write.csv(SH_groups,"SH_groupsB.csv")
#     
# plot(emmeans(SH_Community, ~Year|~Community), comparisons = TRUE) +
#       coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_wrap(~Community)
# 
# SH_Com <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(Community))) +  ggthe_bw +
#   stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1, 
#                    position=position_dodge(width=0.95))+
#   stat_summary(fun.data = "mean_cl_boot",geom = "point", 
#                    position=position_dodge(width=0.95))+
#   stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.95)) +
#   
#   scale_colour_manual(values=c("blue","purple", "red"),
#                         guide = guide_legend(title.position = "top", nrow = 1),
#                         name = "Symbionts detectec in 2014",
#                         labels=c("Cladocopium", "Mixed", "Durusdinium"))+
#   theme(legend.position=c(0.5,0.1),
#         strip.background =element_rect(fill="white")) +
#   scale_y_continuous(limits = c(-3, -0.8),
#             name="\n Total symbiont to host cell ratio (S/H)", 
#             labels = (math_format(10^.x)),
#             expand = c(0.0, 0)) + 
#   scale_x_discrete(name="Year")
# SH_Com

Figure 3: ITS2

SH_Communities_I<-lmerTest::lmer(logTot.SH ~ Year * ITS2_14 * Community + (1|Colony), data=data.Uva)
    summary(SH_Communities_I)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * ITS2_14 * Community + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 80.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2665 -0.4690 -0.1486  0.3796  3.3861 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.02122  0.1457  
##  Residual             0.11122  0.3335  
## Number of obs: 87, groups:  Colony, 29
## 
## Fixed effects:
##                                  Estimate Std. Error       df t value
## (Intercept)                      -1.66014    0.14857 67.55796 -11.174
## Year15                           -0.56741    0.19255 44.59436  -2.947
## Year16                           -0.77804    0.19255 44.59436  -4.041
## ITS2_14C_b                        0.10119    0.25532 70.67197   0.396
## ITS2_14D1                         0.20013    0.51863 70.63905   0.386
## CommunityMixed                   -0.15133    0.47929 70.99075  -0.316
## CommunityDurusdinium-only        -0.04237    0.50922 68.89495  -0.083
## Year15:ITS2_14C_b                 0.67029    0.30501 51.27088   2.198
## Year16:ITS2_14C_b                 0.44786    0.44197 54.72157   1.013
## Year15:ITS2_14D1                  0.72446    0.47149 54.95813   1.537
## Year16:ITS2_14D1                  0.24593    0.56263 57.15914   0.437
## Year15:CommunityMixed             0.20111    0.37100 59.52765   0.542
## Year16:CommunityMixed             0.53106    0.44747 59.20686   1.187
## Year15:CommunityDurusdinium-only -0.08959    0.45989 56.55615  -0.195
## Year16:CommunityDurusdinium-only  0.42021    0.55196 58.74860   0.761
## ITS2_14C_b:CommunityMixed        -0.02664    0.39946 69.55482  -0.067
##                                  Pr(>|t|)    
## (Intercept)                       < 2e-16 ***
## Year15                           0.005090 ** 
## Year16                           0.000208 ***
## ITS2_14C_b                       0.693066    
## ITS2_14D1                        0.700746    
## CommunityMixed                   0.753132    
## CommunityDurusdinium-only        0.933928    
## Year15:ITS2_14C_b                0.032520 *  
## Year16:ITS2_14C_b                0.315370    
## Year15:ITS2_14D1                 0.130144    
## Year16:ITS2_14D1                 0.663686    
## Year15:CommunityMixed            0.589788    
## Year16:CommunityMixed            0.240040    
## Year15:CommunityDurusdinium-only 0.846244    
## Year16:CommunityDurusdinium-only 0.449526    
## ITS2_14C_b:CommunityMixed        0.947010    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 11 columns / coefficients
    #step(SH_Communities_I)
    
    SH_CI.emmc<-emmeans(SH_Communities_I, ~Year*ITS2_14*Community)
    SH_CI_groups<-cld(SH_CI.emmc)
    SH_CI_groups<-as.data.frame(SH_CI_groups[complete.cases(SH_CI_groups),])
    SH_CI_groups<-SH_CI_groups[order(SH_CI_groups$Year,SH_CI_groups$ITS2_14, SH_CI_groups$Community),]
    SH_CI_groups
    #write.csv(SH_groups,"SH_Communities_multicomp.csv")

SH_Comties_I <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(Community))) + ggthe_bw+
    
      facet_grid(~ITS2_14, labeller = labeller(ITS2_14=c(
                      `C_a` = "C_a (2014)",
                      `C_b` = "C_b (2014)", 
                      `D1` = "D1 (2014)"))) +
        
      stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
      stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
      stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
      
      scale_colour_manual(values=c("blue3","purple", "red3"),
                        guide = guide_legend(title.position = "top", nrow = 3),
                        name = "Detected symbionts",
                        labels=c("Cladocopium", "Mixed", "Durusdinium"))+
      theme(legend.position=c(0.82,0.2),
            legend.box.background = element_rect(),
            strip.background =element_rect(fill="white"))+
   
      scale_y_continuous(limits = c(-3, -0.8),
            name="\n Total symbiont to host cell ratio (S/H)", 
            labels = (math_format(10^.x)),
            expand = c(0.0, 0)) + 
      scale_x_continuous(name="", 
                         breaks = c(2014, 2015, 2016),
                         labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) 
SH_Comties_I

#ggsave(file="Figure_3.svg", plot=SH_Comties, width=6, height=4)

8b. Changes in Cladocopium abundance

CH_Com <- ggplot(data.Uva, aes (Year.2, logC.SH, colour=factor(Community))) + ggthe_bw+ ggtitle("a.")+
    
      stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
      stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
      stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
      
      scale_colour_manual(values=c("blue3","purple", "red3"),
                        guide = guide_legend(title.position = "top", nrow = 3),
                        name = "Detected symbionts",
                        labels=c("Cladocopium", "Mixed", "Durusdinium"))+
      theme(legend.position=c(0.82,0.8),
            legend.box.background = element_rect(),
            strip.background =element_rect(fill="white"))+
   
      scale_y_continuous(limits = c(-3, -0.7),
            name="\n Cladocopium to host cell ratio (C/H)", 
            labels = (math_format(10^.x)),
            expand = c(0.0, 0)) + 
      scale_x_continuous(name="", 
                         breaks = c(2014, 2015, 2016),
                         labels = c(" ", " ", " ")) 
#CH_Com

CH_Com2ITS2<-CH_Com + facet_grid(~ITS2_14)
# CH_Com2ITS2

CH_I<-lmer(logC.SH ~ Year* Community* ITS2_14 + (1|Colony), data=data.Uva)
    summary(CH_I)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logC.SH ~ Year * Community * ITS2_14 + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 107.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.72727 -0.28740 -0.04472  0.44178  1.54498 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.1101   0.3318  
##  Residual             0.3945   0.6281  
## Number of obs: 55, groups:  Colony, 21
## 
## Fixed effects:
##                       Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)           -1.66014    0.29000 39.38253  -5.725 1.21e-06 ***
## Year15                -0.56741    0.36265 23.92752  -1.565   0.1308    
## Year16                -0.77804    0.36265 23.92752  -2.145   0.0423 *  
## CommunityMixed        -0.38598    0.48624 40.68407  -0.794   0.4319    
## ITS2_14C_b             0.23717    0.49523 42.78435   0.479   0.6344    
## ITS2_14D1             -0.91902    0.63611 42.99226  -1.445   0.1558    
## Year15:CommunityMixed  0.31035    0.70821 33.47664   0.438   0.6640    
## Year16:CommunityMixed -0.34608    0.85164 32.26324  -0.406   0.6872    
## Year15:ITS2_14C_b      0.55621    0.57741 27.86114   0.963   0.3437    
## Year16:ITS2_14C_b      0.36653    0.83852 29.64957   0.437   0.6652    
## Year15:ITS2_14D1       0.04458    0.89702 30.76732   0.050   0.9607    
## Year16:ITS2_14D1       1.23468    1.06989 31.29787   1.154   0.2572    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##              (Intr) Year15 Year16 CmmntM ITS2_14C ITS2_14D Y15:CM Y16:CM
## Year15       -0.625                                                     
## Year16       -0.625  0.500                                              
## CommuntyMxd   0.000  0.000  0.000                                       
## ITS2_14C_b   -0.586  0.366  0.366 -0.655                                
## ITS2_14D1    -0.456  0.285  0.285 -0.764  0.767                         
## Yr15:CmmntM   0.000  0.000  0.000 -0.628  0.411    0.480                
## Yr16:CmmntM   0.000  0.000  0.000 -0.501  0.328    0.383    0.330       
## Y15:ITS2_14C  0.393 -0.628 -0.314  0.546 -0.740   -0.596   -0.508 -0.278
## Y16:ITS2_14C  0.270 -0.216 -0.432  0.315 -0.470   -0.364   -0.219 -0.810
## Y15:ITS2_14D  0.253 -0.404 -0.202  0.496 -0.473   -0.610   -0.790 -0.261
## Y16:ITS2_14D  0.212 -0.169 -0.339  0.399 -0.385   -0.498   -0.263 -0.796
##              Y15:ITS2_14C Y16:ITS2_14C Y15:ITS2_14D
## Year15                                             
## Year16                                             
## CommuntyMxd                                        
## ITS2_14C_b                                         
## ITS2_14D1                                          
## Yr15:CmmntM                                        
## Yr16:CmmntM                                        
## Y15:ITS2_14C                                       
## Y16:ITS2_14C  0.404                                
## Y15:ITS2_14D  0.655        0.260                   
## Y16:ITS2_14D  0.327        0.791        0.352      
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 6 columns / coefficients
    #step(CH_I)

    CH_I.emmc<-emmeans(CH_I, ~Year*Community* ITS2_14)
    CH_I_groups<-cld(CH_I.emmc)
    CH_I_groups<-as.data.frame(CH_I_groups[complete.cases(CH_I_groups),])
    CH_I_groups<-CH_I_groups[order(CH_I_groups$Year,CH_I_groups$ITS2_14, CH_I_groups$Community),]
    CH_I_groups
    #write.csv(CH_I_groups,"CH_I_groups.csv")

9b. Changes in Durusdinium abundance

DH_ComI <- ggplot(data.Uva, aes (Year.2, logD.SH, colour=factor(Community))) + 
      ggthe_bw+ ggtitle("b.")+
    
      stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
      stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
      stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
      
      scale_colour_manual(values=c("purple", "red3"),
                        guide = guide_legend(title.position = "top", nrow = 3),
                        name = "Detected symbionts",
                        labels=c("Mixed", "Durusdinium"))+
      theme(legend.position=c(0.82,0.2),
            legend.box.background = element_rect(),
            strip.background =element_rect(fill="white"))+
   
      scale_y_continuous(limits = c(-6, -0.7),
            name="\n Durusdinium to host cell ratio (D/H)", 
            labels = (math_format(10^.x)),
            expand = c(0.0, 0)) + 
      scale_x_continuous(name="", 
                         breaks = c(2014, 2015, 2016),
                         labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016"))

DH_Com2ITS2<-DH_ComI + facet_grid(~ITS2_14)
DH_Com2ITS2

DH_LM_ITS2<-lmer(logD.SH ~ Year* Community* ITS2_14 + (1|Colony), data=data.Uva)
    summary(DH_LM_ITS2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logD.SH ~ Year * Community * ITS2_14 + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 66.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0953 -0.5756  0.0079  0.3844  3.6310 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.0000   0.0000  
##  Residual             0.1687   0.4107  
## Number of obs: 58, groups:  Colony, 22
## 
## Fixed effects:
##                                     Estimate Std. Error       df t value
## (Intercept)                         -4.86573    0.16766 48.00000 -29.021
## Year15                               0.07049    0.33532 48.00000   0.210
## Year16                               3.21459    0.23711 48.00000  13.557
## CommunityDurusdinium-only            0.27806    0.50992 48.00000   0.545
## ITS2_14D1                            3.15614    0.23711 48.00000  13.311
## Year15:CommunityDurusdinium-only    -0.35506    0.32897 48.00000  -1.079
## Year16:CommunityDurusdinium-only    -0.22094    0.38416 48.00000  -0.575
## Year15:ITS2_14D1                     0.36688    0.42745 48.00000   0.858
## Year16:ITS2_14D1                    -3.09697    0.41068 48.00000  -7.541
## CommunityDurusdinium-only:ITS2_14D1 -0.08557    0.45916 48.00000  -0.186
##                                     Pr(>|t|)    
## (Intercept)                          < 2e-16 ***
## Year15                                 0.834    
## Year16                               < 2e-16 ***
## CommunityDurusdinium-only              0.588    
## ITS2_14D1                            < 2e-16 ***
## Year15:CommunityDurusdinium-only       0.286    
## Year16:CommunityDurusdinium-only       0.568    
## Year15:ITS2_14D1                       0.395    
## Year16:ITS2_14D1                    1.09e-09 ***
## CommunityDurusdinium-only:ITS2_14D1    0.853    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Year15 Year16 CmmnD- ITS2_1 Y15:CD Y16:CD Y15:IT Y16:IT
## Year15      -0.500                                                        
## Year16      -0.707  0.354                                                 
## CmmntyDrsd-  0.000  0.000 -0.232                                          
## ITS2_14D1   -0.707  0.354  0.500 -0.232                                   
## Yr15:CmmnD-  0.000  0.000  0.000 -0.293  0.360                            
## Yr16:CmmnD-  0.000  0.000  0.000 -0.753  0.309  0.389                     
## Y15:ITS2_14  0.392 -0.784 -0.277  0.129 -0.555 -0.500 -0.171              
## Y16:ITS2_14  0.408 -0.204 -0.577  0.671 -0.577 -0.208 -0.713  0.320       
## CD-:ITS2_14  0.000  0.000  0.258 -0.900  0.000  0.000  0.558  0.000 -0.596
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
## convergence code: 0
## boundary (singular) fit: see ?isSingular
    #step(DH_LM_ITS2)

library(emmeans)
    DH_I.emmc<-emmeans(DH_LM_ITS2, ~Year*Community* ITS2_14)
    DH_groups_I<-cld(DH_I.emmc)
    DH_groups_I<-as.data.frame(DH_groups_I[complete.cases(DH_groups_I),])
    DH_groups_I<-DH_groups_I[order(DH_groups_I$Year,DH_groups_I$ITS2_14, DH_groups_I$Community),]
    DH_groups_I
    #write.csv(DH_groups_I,"DH_I_groups.csv")

Figure 4: qPCR

Figure_4b<-grid.arrange(CH_Com2ITS2, DH_Com2ITS2,nrow=2)

# ggsave(file="Outputs/Figure_4b.svg", plot=Figure_4, width=7, height=4)

10. Supplmenmtary C/H D/H dynamics

Select data

log_Log.data<-data.Uva[,c("Colony","Year.2", "C.SH", "D.SH", 
                          "D14Community", "Community","ITS2_14")]

log_Log.data<-log_Log.data %>%
  filter(Year.2 != "2015")

log_Log.data$logC<-log10(log_Log.data$C.SH)
log_Log.data$logD<-log10(log_Log.data$D.SH)
log_Log.data$logC[which(log_Log.data$C.SH==0)] <- -6
log_Log.data$logD[which(log_Log.data$D.SH==0)] <- -6

Plot Black and White

Community_Changes2 <- ggplot (log_Log.data, aes(logC, logD, gropup=Colony)) +

  geom_abline(intercept=0, slope=1, linetype=1, color="black")+
  geom_abline(intercept=-1, slope=1, linetype=2, color="black")+
  geom_abline(intercept=-2, slope=1, linetype=2, color="darkgray")+
  geom_abline(intercept=-3, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=-4, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=-5, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=1, slope=1, linetype=2, color="black")+
  geom_abline(intercept=2, slope=1, linetype=2, color="darkgray")+
  geom_abline(intercept=3, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=4, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=5, slope=1, linetype=2, color="gray")+
  
  ggthe_bw + 
  scale_x_continuous(limits = c(-6,-0.1), "log10 (C/H cell ratio)")+
  scale_y_continuous(limits = c(-6,-0.1), "log10 (D/H cell ratio)") +
  
  geom_point(alpha=0.5, size=4) +
  geom_path(size = 1,arrow =
              arrow(type = "open", angle = 30, length = unit(0.1, "inches")))+
  
  annotate("text", x = c(seq(-5.1, -1.1, by=1)), y = -0.9, label = c(
    "10,0000C : D", "1,000C : 1D", "100C : 1D","10C : 1D", "C : D"), angle = 45)+
  annotate("text", x = -0.9, y = c(seq(-4.7, -1.7, by=1)), label = c(
    "1C : 10,000D", "1C : 1,000D", "1C : 100D",  "1C : 10D"), angle = 45)
  
Community_Changes2

#ggsave(file="FigAndrewsFavoriteGraph_BW.svg", plot=Community_Changes2, width=5, height=5)
Community_Changes3 <- ggplot (log_Log.data, aes(logC, logD, gropup=Colony, colour=D14Community)) +

  geom_abline(intercept=0, slope=1, linetype=1, color="black")+
  geom_abline(intercept=-1, slope=1, linetype=2, color="black")+
  geom_abline(intercept=-2, slope=1, linetype=2, color="darkgray")+
  geom_abline(intercept=-3, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=-4, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=-5, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=1, slope=1, linetype=2, color="black")+
  geom_abline(intercept=2, slope=1, linetype=2, color="darkgray")+
  geom_abline(intercept=3, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=4, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=5, slope=1, linetype=2, color="gray")+
  scale_colour_manual(values=c("blue","purple", "red"),
                        guide = guide_legend(title.position = "top", nrow = 3),
                        name = "Detected symbionts (2014)",
                        labels=c("Cladocopium", "Mixed", "Durusdinium"))+
  
  ggthe_bw + theme(legend.position="bottom")+
  scale_x_continuous(limits = c(-6,-0.1), "log10 (C/H cell ratio)")+
  scale_y_continuous(limits = c(-6,-0.1), "log10 (D/H cell ratio)") +
  
  geom_point(alpha=0.5, size=4) +
  geom_path(size = 1,arrow =
              arrow(type = "open", angle = 30, length = unit(0.1, "inches")))+
  
  annotate("text", x = c(seq(-5.1, -1.1, by=1)), y = -0.9, label = c(
    "10,0000C : D", "1,000C : 1D", "100C : 1D","10C : 1D", "C : D"), angle = 45)+
  annotate("text", x = -0.9, y = c(seq(-4.7, -1.7, by=1)), label = c(
    "1C : 10,000D", "1C : 1,000D", "1C : 100D",  "1C : 10D"), angle = 45)
  
Community_Changes3

# ggsave(file="FigAndrewsFavoriteGraph_qPCR.svg", plot=Community_Changes3, width=6, height=6)
Community_Changes4 <- ggplot (log_Log.data, aes(logC, logD, gropup=Colony, colour=ITS2_14)) +

  geom_abline(intercept=0, slope=1, linetype=1, color="black")+
  geom_abline(intercept=-1, slope=1, linetype=2, color="black")+
  geom_abline(intercept=-2, slope=1, linetype=2, color="darkgray")+
  geom_abline(intercept=-3, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=-4, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=-5, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=1, slope=1, linetype=2, color="black")+
  geom_abline(intercept=2, slope=1, linetype=2, color="darkgray")+
  geom_abline(intercept=3, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=4, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=5, slope=1, linetype=2, color="gray")+
  scale_colour_manual(values=c("blue","purple", "red"),
                        guide = guide_legend(title.position = "top", nrow = 3),
                        name = "Main ITS2 types (2014)",
                        labels=c("C_a", "C_b", "D1"))+
  
  ggthe_bw + theme(legend.position="bottom", 
                    legend.box = "horizontal")+
  scale_x_continuous(limits = c(-6,-0.1), "log10 (C/H cell ratio)")+
  scale_y_continuous(limits = c(-6,-0.1), "log10 (D/H cell ratio)") +
  
  geom_point(alpha=0.5, size=4) +
  geom_path(size = 1,arrow =
              arrow(type = "open", angle = 30, length = unit(0.1, "inches")))+
  
  annotate("text", x = c(seq(-5.1, -1.1, by=1)), y = -0.9, label = c(
    "10,0000C : D", "1,000C : 1D", "100C : 1D","10C : 1D", "C : D"), angle = 45)+
  annotate("text", x = -0.9, y = c(seq(-4.7, -1.7, by=1)), label = c(
    "1C : 10,000D", "1C : 1,000D", "1C : 100D",  "1C : 10D"), angle = 45)
  
Community_Changes4

#ggsave(file="FigAndrewsFavoriteGraph_ITS2.svg", plot=Community_Changes4, width=6, height=6)